  //__syncthreads();
 
  /*if(threadIdx.x==0 && threadIdx.y==0)
  {
     //Boundary case Y = 0 
     for(int x=0; x<BLOCKDIM_X; x++)
     {
         Eprev_shared[0][x] = Eprev[index - width + x];
     } 
     //Boundary case Y = BLOCKDIM_Y-1
     for(int x=0; x<BLOCKDIM_X; x++)
     {
         Eprev_shared[BLOCKDIM_Y+1][x] = Eprev[index + width*BLOCKDIM_X + x];
     } 
  }
  if(threadIdx.x==BLOCKDIM_X-1 && threadIdx.y==BLOCKDIM_Y-1)
  {
     //Boundary case X = 0 
     for(int y=0; y<BLOCKDIM_Y; y++)
     {
         //Eprev_shared[y][BLOCKDIM_Y+1] = Eprev[index + BLOCKDIM_Y + width*y + 1];
     } 
     //Boundary case X = BLOCKDIM_X-1
     for(int y=BLOCKDIM_Y; y<=0; y--)
     {
         Eprev_shared[BLOCKDIM_Y+1][y] = Eprev[index - width*y + 1];
     }

  }*/

  //load data elements for thread in on the 
  //boundary of the shared memory tile
  //load data items on boundary x=0  

__global__ void boundaries(_DOUBLE_ *Eprev, const int n, const int m)
{
  int width = n + 3;

  int idx= threadIdx.x + 1;
  int idy= threadIdx.y + 1;

  int gidx = idx + blockDim.x * blockIdx.x ;
  int gidy = idy + blockDim.y * blockIdx.y ;

  if(gidx > n+1 || gidy > m+1)
   return;
 
  int index = gidy * width + gidx ; 

  __syncthreads();
  if (gidx == 2 )
  {
    Eprev[index-2] = Eprev[index];
    Eprev[(index+n-2)+2] = Eprev[(index+n-2)];
  }

  __syncthreads();
  if(gidy == 2 )
  {
     Eprev[index - 2*width] = Eprev[index];
     Eprev[index + (m-2)*width + 2*width ] = Eprev[index + (m-2)*width];
  }
}
// Collect the boundary data
__global__ void boundariesY(_DOUBLE_ *Eprev, const int n, const int m)
{
  int width = n + 3;

  int idx= threadIdx.x +1;
  int idy= threadIdx.y +1;

  int gidx = idx + blockDim.x * blockIdx.x ;
  int gidy = idy + blockDim.y * blockIdx.y ;

if(gidx > n+1 || gidy > m+1)
  return;
 
  int index = gidy * width + gidx ; 
   __syncthreads(); 
  if(gidy == 2 )
    Eprev[index - 2*width] = Eprev[index] ;
   __syncthreads(); 
  if ( gidy == m)
    Eprev[index + 2*width ] = Eprev[index];

}
__global__ void boundariesX(_DOUBLE_ *Eprev, const int n, const int m)
{
  int width = n + 3;

  int idx= threadIdx.x +1;
  int idy= threadIdx.y +1;

  int gidx = idx + blockDim.x * blockIdx.x ;
  int gidy = idy + blockDim.y * blockIdx.y ;

if(gidx > n+1 || gidy > m+1)
  return;
 
  int index = gidy * width + gidx ; 

   __syncthreads(); 
  if (gidx == 2 )
    Eprev[index-2] = Eprev[index];
   __syncthreads(); 
  if(gidx == n)
    Eprev[index+2] = Eprev[index];

}


  //__syncthreads();
  /*if (gidy >= 1 && gidy <= m+1) {
	  if (gidx >= 1 && gidx <= n+1) {
	    E[index] = Eprev[index]+alpha*(Eprev[index+1]+Eprev[index-1]
					    -4*Eprev[index]+Eprev[index + n + 3]+Eprev[index - n - 3]);
       E[index] += -dt*(kk*E[index]*(E[index]-a)*(E[index]-1)+E[index]*R[index]);
       R[index] += dt*(epsilon+M1* R[index]/( E[index]+M2))*(-R[index]-kk*E[index]*(E[index]-b-1));
          }
  }*/
  /*__syncthreads();
  if (gidy >= 1 && gidy <= m+1) {
	  if (gidx >= 1 && gidx <= n+1) {
	    E[index] += -dt*(kk*E[index]*(E[index]-a)*(E[index]-1)+E[index]*R[index]);
          }
  }
  __syncthreads();
  if (gidy >= 1 && gidy <= m+1) {
	  if (gidx >= 1 && gidx <= n+1) {

	    R[index] += dt*(epsilon+M1* R[index]/( E[index]+M2))*(-R[index]-kk*E[index]*(E[index]-b-1));
	  }
    }*/
